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Abstract. Two effectively one-dimensional parallel coupled Bose-Einstein conden- 
sates in the presence of external potentials are studied. The system is modelled by 
linearly coupled Gross-Pitacvskii equations. In particular, grey-soliton-like solutions 
representing analogues of superconducting Josephson fluxons as well as coupled dark 
solitons are discussed. Theoretical approximations based on variational formulations 
are derived. It is found that the presence of a magnetic trap can destabilize the fluxon 
analogues. However, stabilization is possible by controlling the effective linear coupling 
between the condensates. 
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1. Introduction 

The concept of electron tunnelling between two superconductors separated by a thin 
insulating barrier predicted by Josephson [1] has been extended relatively recently to 
tunnelling of Bose-Einstein condensates (BECs) across a potential barrier by Smerzi et 
al. [21 131 H]. Such tunnelling has been observed experimentally where a single [5l E] 
and an array [7] of short Bose- Josephson junctions (BJJs) were realized. The dynamics 
of the phase difference between the wavefunctions of the condensates [21 [3l HI El El [ID] 
resembles that of point-like Josephson junctions 

Recently a proposal for the realization of a long BJJ has been presented by Kaurov 
and Kuklov [121 IIS]- Similarly to superconducting long Josephson junctions, one may 
also look for an analogue of Josephson fluxons [14] in this case. It was shown in [121 113] 
that fluxon analogues are given by coupled dark-soliton-like solutions, as the relative 
phase of the solutions has a kink shape with the topological phase difference equal to 
27r. Moreover, it was emphasized that fluxon analogues (FAs) can be spontaneously 
formed from coupled dark solitons due to the presence of a critical coupling at which 
the two solitonic structures exchange their stability. The idea of FAs in tunnel-coupled 
BECs is then extended to rotational FAs in the ground state of rotating annular BECs 
confined in double-ring traps |15j . 

In this report, we consider the existence and the stability of FAs in two coupled 
cigar-shaped condensates in the presence of a magnetic trap along the elongated 
direction modelled by the normalized coupled Gross-Pitaevskii equations [121 US] 



where ifjj, j = 1, 2, is the bosonic field, and t and x are the time and axial coordinate, 
respectively. Here, we assume that the parallel quasi one-dimensional BECs are linked 
effectively by a weak coupling k. Note that herein fc > 0. The case k < can be 
obtained accordingly as there is a symmetry transformation k — )■ —k and ijjj — )■ iipj. 
The parameters /i and u are respectively the nonlinearity coefficient representing the 
atomic scattering length and the chemical potential. Here, we consider the defocusing 
nonlinearity /i < 0, see, e.g., [Ml [HI [181 IlH]) for the focusing counterpart /i > 0. V is 
the magnetic trap with strength Q, i.e. 



Extending the work [121 which was without the magnetic trap, i.e. f2 = 0, we are 
particularly interested in the effects of > to the localized excitations. 

When f2 = 0, writing ipj = \ipj\ex]){iifj), it was shown that the relative phase 
= ±(<y52 — ifi) will satisfy a modified sine-Gordon equation [12]. A fluxon analogue of 
([1]) in that case is given by the solution xl^i = ip2 = with 



(1) 




(2) 




(3) 
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where the asterisk denotes complex conjugation. The sohton can be regarded as an 
analogue of Josephson fiuxons [121 US] as the phase difference between the phases of ipi 
and ip2 forms a spatial kink connecting = and = ±27r. In the following, solution 
(|3]) (and its continuations) will be referred to as FAs. From the expression, it is clear 
that an FA exists only for < k < uj/3. For k = uj/3, the solution in (jS]) transforms 
into dark solitons [121 US] 

■01^2 = ±^ - tanh(-\/u; + fcx). (4) 

which exists for k > — w. Thus solutions in (|3|) and (|4|) coincide at k = u/3 and hence 
k = u/3 is the bifurcation point along the family (jl]). In the following, we denote this 
critical coupling as kce- 

When the two condensates are uncoupled, i.e. k = 0, the dynamics of a dark 
soliton in BECs with magnetic trap has been considered before theoretically [20| [22] 
(see also [23] and references therein) and experimentally [241 EEl [271 ESI [28] . Interesting 
phenomena on the collective behavior of a quantum degenerate bosonic gas, such as 
soliton oscillations [271 ISSl [21] and frequency shifts due to soliton collisions [28] were 
observed. A theoretical analysis based on variational formulation was developed in 
[291 [30] that is in good agreement with numerics as well as with experiments (see, e.g., 
[311 [32]). A similar variational method will be derived here to explain the dynamics of 
FAs. 

The present report is outlined as follows. In Section 2 we will derive a relation 
between the velocity of the travelling FAs and the critical value kce of the coupling 
constant at which the FAs solution changes into a coupled dark soliton. In Section 3, 
a variational formulation for the FAs is considered to study their dynamical behaviours 
analytically. We will solve the governing equations numerically in Section 4. Comparing 
the numerical and the analytical results, we show good agreement between them. In 
the last section, we conclude the work and present open problems for future work. 



2. The velocity dependence of the critical coupling k^e 

When FAs do not travel, it has been discussed above that kce = u}/3. When the solitons 
move with velocity v, it is natural to expect that kce = kce{v). It is because dark solitons 
only exist for f < 1 [29], while at kce FAs become dark solitons. In the following, we 
are interested in the expression of the critical coupling. We will determine the velocity 
dependence of kce analytically. 

For simplicity, first we scale the governing equation ([1]), such that the nonlinearity 
coefficient and the chemical potential become n = u = 1. We also scale the wavefunction 
ijjj in ([1]) by = ipj/ \/l + k, such that the equation becomes 

2 ' ^' 1 + A; 1 + k' ^ ^ 



where t = (1 + k)t and x = \/l + k x. 
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Defining ^ 
written as 



X 



vt, Eq. (E]) with V = in a moving coordinate frame can be 



1 



3-i 



k 



0. 



(6) 
(7) 



Tfie equation above lias an explicit solution for the travelling dark soliton 

ipds = Atanh(5^) + iv, 

with A = B and + = 1. 

Let us perturb ipds such that ipj = i^ds + {~^)^"^'^f{0^ where e ^ 1. Substituting 
these values of tpi and 1^2 in (El) and retaining linear terms in e only, we obtain 

1 d'fiO 



+ 



2 de 
2 J. 1,2 



2{AHanh\B^) + v'^) + k 



tanr(5{) 



no 



2Ai;tanh(E0/*(f) 



dl J 



0, 



where we have substituted v = \/l + kv. Here, v is the velocity measured in the 'original' 
time t, while v is in the scaled time t. 
Choosing 

/(O = p sech(Sf) tanh(EO + « q sech(SO, (8) 

where p and q are real, and substitute it in the last equation above, we obtain two 
equations containing p and q. Eliminating p and q in the resulting equations, we obtain 
the following relation between k and v 



k 



1 

21 



4 
21 



V7t;4 - 7?;2 + 4. 



(9) 



The coupling constant k here is the critical coupling kce, which is clearly a function of 
V. Note from the above expression that fcce — >■ | as u — )■ and fcce — >■ as w — )■ 1. 



3. Variational formulations of the FAs 

The Lagrangian formalism for the FAs will be separated into two cases, i.e. when the 
coupling constant k is close to the critical coupling kce at which the FAs are close to 
coupled dark solitons and when it is close to the uncoupled limit k ^ 0. The two cases 
determine the ansatz that will be taken to approximate the solutions. 



3.1. The case k ~ kce md V ^ 

When there is a harmonic potential, the dynamics of the solutions will be determined 
perturbatively. We will particularly follow the method of |29l [SI |32], which was 
developed for dark solitons, to our case. A similar calculation will be performed to 
discuss the dynamics and stability of FAs in the limit of k close to kce as in that case 
the real part, i.e. the dark soliton component, is dominant and Ci, C2 ~ 0. 
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Figure 1. (a) Numerically obtained coupled dark solitons for a; ~ 1. fi ~ 1, and 
k = 0.1. The black curves represent the real parts while the horizontal red lines are 
the imaginary parts of ipj , respectively, (b) The eigenvalues A of the solutions in (a) in 
the complex plane showing the instability of the solutions, (c) The critical eigenvalue 
of coupled dark solitons as a function of k. (d) An evolution of the unstable coupled 
dark solitons in (a). Shown is 1-0^1. A grcy-soliton-likc structure in the centre at the 
end of the computation is an FA as the spatial profile of the phase difference between 
^pl and ip2 form a 27r-kink shape (not shown here). 



Due to the presence of the magnetic trap, which is assumed to be slowly varying, 
i.e. fi^ ^ 1, we take the ansatz 

= ^TF$i, (10) 

where "^tf is the Thomas- Fermi cloud approximately given by 

^TF = Vmax{(l -1^/(1 + /c)),0}. (11) 

and 

$j = Atanh (z) + i {Cj sech (z) C) , j = 1, 2, (12) 

where z = B [x — xq). The parameters A, B, C, xq are in general functions of t. Here, 
+ v"^ = 1. When Cj = 0, one can recognise that the above function is the usual 
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Figure 2. (a) A numerically obtained FAs for = 1, /i = 1, and k = 0.1. The black 
curves represent the real parts while the red curves are the imaginary parts of "01 and 
^p2, respectively, (b) The eigenvalue distribution of the FAs in (a) in the complex plane 
showing the stability of the solution. 



ansatz to describe a dark soliton travelling with the velocity v, where in that case the 
parameters A, B, v and Xq can be conveniently written as A = B = cos(f, v = simp, 
and Xo(i) = J^siinp^t') dt' [3T] . 

Note the similarity of ffTOj) with the ansatz = '^xF'^^ds ISH [32] used to study the 
dynamics of dark solitons in a harmonic potential. Substituting the ansatz ffTOj) into ([5]) 
yields the equation for $j, i.e. 



Jt 2 ■' 



3-J 



1 + k 



(13) 



where 



1 



1, 



[l + k) 

Here, for R we only retain terms linear in V and V^, but neglecting terms linear in Vxx- 
When there is no magnetic trap, (|T3l) with V = Q can be derived from the 
Lagrangian 

2 



1*7- 



k) 



1 + k 



dx. 



(14) 



It is not clear how to evaluate the first term of the integral ( 1141) due to the present of 



nonzero in the denominator of (1 — 1/ 



Because of that, in the following we 



assume Cj to be small, i.e. we consider the case of /c ~ kce, and take the series expansion 
with respect to Cj about Cj = upto 0(C|). Performing the integration will yield the 
effective Lagrangian 



- 4Av + 4 tan" 



-Hi)-7rAC, 



A\A' + B' 
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B{l + k) 



1 

35 



_ 16^2)^4 + 2^6 



(3A; + 2)^4 - 2A;CiC2 + 0{C: 



^3 I-' 



(15) 



where C^^^ = (Ci)'^ + (C2)^ k 



1, . . . , 4. One can check that when the sohton is not 
moving, i.e. Xq^ = v = and A = 1, the Euler-Lagrange equations for the remaining 
parameters B, Ci, and C2 can be solved analytically to yield 

2Vk 



5(0) 



VTTk 



CI 



(0) 



VI -3k 



or 



" VTTk 
5(0) ^ 1 c*!") 



(16) 

0, which is the dark 



which is nothing else, but the solution ^ ^ ^' 
soliton (jl]). It is important to note here that despite the series expansion with respect 
to Cj that we took prior to integrating ( fT4|) . our result ( TT6l) corresponds to an exact 
solution for any value of < < kce- It is not yet clear to us why this is the case. 

To determine the influence of the magnetic trap, i.e. 7^ 0, to the FAs, we will 
follow the idea of [22] and treat the right hand side of (IT^ as perturbations. Due to 
the presence of the perturbations the Euler-Lagrange equation for the variable aj, with 



cp, B, and Cj, becomes 



dC. 



dC 



2 Re 




Following [211 [22], we need two active variables only. Because of that, we will assume 
that adiabatically B and Cj are independent of time and are given by (fT6l) . Hence, we 
obtain the following set of equations 

dA VL^xqv 



di l2Vk{k + l)-2 
dxf) V 
di ~ 57QAkl{k + l)l 

+ 576k^{A^ + 2) + 



A'{k 



Ilk 



(17) 



192k\A^ + 14) + 576^(^2 
n^e\^{4Sxo^ + 7i^){A^ + 3) 



+ 12(^2 + 9)| + 2A:fi2 1 24x0^(^2 - 1) + 12(^2 + 4) 
+ 71^(^2 + 1)} + n\A^ - l)(7r2 + 12) + 192k{A^ - 4) 

From the two equations above, we obtain (see, e.g., [32] ) 
d'xo (1 - 5k)Q' 

-Xq + U[\l ), 



:i8) 



(19) 



dt^ 1 + k 

where we have assumed that A ^ 1 and a;o ~ 0. Note that the derivative in f|T9l) 
is with respect to the original time variable. This shows that the solution is stable for 
values of k > ^. When /c = 1/3, we recover the oscillation frequency of dark solitons in a 
harmonic trap [21] (see also [31], [221 [21] ) . We will check the validity of our approximation 
by comparing it with numerical results. 
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Figure 3. (a) A numerically obtained FAs for v ~ 0.2 and k ~ 0.1. (b) The critical 
coupling constant fc^e as a function of the velocity v for the existence of travelling FAs. 
Filled circles are numerical data and the solid line is the function in ({9]). FAs only exist 
on the left of the solid curve. 



3.2. The case k^O andV ^0 

When k is close to zero, the imaginary part of the FAs is dominant over the real part. 
Due to the presence of a magnetic trap, the imaginary part will become the Thomas- 
Fermi ground state in the limit — )■ 0. Motivated by, e.g., [33] where a Gaussian ansatz, 
which is the ground state solution of the linear Schrodinger equation ([1]) and ([2]) with 
/i = /c = 0, was shown to be able to approximate the ground state of the nonlinear 
equation, in the following we will take a similar ansatz for our problem. For the present 
case, we will consider the governing equation ([1]) to be derived from the Lagrangian 

+ 2 {cu - V) + 2ktpjij;_j dx. (20) 

Note that the main difference between (!20|) and (IMI) is the factor ^1 — prp^ in the first 
term of the integrand. As for the ansatz, we take 

= {AjX + iCj) e-^""'. (21) 

One advantage of such a Gaussian ansatz is that the Lagrangian can be evaluated rather 
straightforwardly. This advantage is expected to be particularly useful in studying 
bound states of FAs. 

Substituting the ansatz into the Lagrangian fl20l) yields the effective Lagrangian 



4^2^3(1252 + 3fi2 - StuB) 

256i?2 L 

+ 16V2BC4{4:B^ + - 8luB) + GAB^Cq 

+ 16B{Aq - 4:V2kA7) + 3^4] , (22) 
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from which we obtain the Euler- Lagrange equations 



sV2BkQ 



3-i + Cj 



0, 



0, (23) 
(24) 



A/ + 72^2 + aV2B{B -2u + V2C/ 

15^4 + 165(3^6 + 45^6) + I6V2BC4 {Sn^ - 452 - 8Buj) 

+ 12^2^3 {5n^ + 452 - 8Buj) - uV2kB {ABC1C2 + 3A1A2) = 0, (25) 

where = Al + Al A^ = A\ + Al A^ = A^C^ + A^ = AjCf + A^Cl 

A'j = A1A2 + 45C1C2. Solving the system of algebraic equations above yields an 
approximation to the FAs. The validity of the ansatz will be checked by comparing the 
results presented here with numerical results. The comparison is presented in section 
14.21 As for the stability, using a Gaussian ansatz, one will need to include, e.g., chirp 
variables, which we leave for future investigation. 



4. Numerical computations 

To study the existence of static FAs in the time-independent framework of ([T]), we use a 
Newton- Raphson continuation method with an initial ansatz ([3]) or (jl]). The spatial first 
and second order derivatives (the former being in the governing equation in a moving 
coordinate frame) are approximated using central finite difference with three-point or 
five-point stencils. At the computational boundaries, we use the Neumann boundary 
conditions. Numerical linear stability analysis of a solution t/jj^\x) is then performed 
by looking for perturbed solutions of the form 

= V^f (^) + eK-(^) e''' + b*{x) J = 1, 2. 

Substituting the ansatz into the governing equation ([1]) and keeping the linear terms 
in e, one will obtain a linear eigenvalue problem for the stability of 'ipf^- The ensuing 
eigenvalue problem is then discretized using a similar finite difference scheme as above 
and solved numerically for the eigenfrequency A and corresponding eigenfunctions aj 
and bj. It is then clear that ip^^\x) is a stable solution if the imaginary parts of all the 
eigenvalues vanish, i.e. Im(A) = 0. 

The evolutions in time of the solitons when they are unstable are investigated 
through integrating the time-dependent governing equations using a Runge-Kutta 
method of order four. 



4.1. Coupled BECs without a trap: = 

As pointed out in [121 [13] for k small enough, coupled dark solitons (j4]) are unstable. 
Shown in Figs. 1(a) and 1(b) are respectively numerically obtained coupled solitons and 
their eigenvalue structure in the complex plane for k = 0.1. Because there is at least 
one pair of eigenvalues with nonzero imaginary parts, we conclude that the solitons are 
unstable. We have calculated the stability of dark solitons for different values of coupling 
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constant k, where we obtained that dark sohton solutions are stable for k > 1/3 as shown 
in Fig. 1(c) In the instability region of the coupled dark solitons, one can obtain stable 



FAs [111 [13]. Shown in Fig. 1(d) is the typical time-dynamics of unstable coupled dark 
solitons, where an FA is obtained spontaneously. Nonetheless, it is important to note 
that it is not always the case. When the coupling constant is too small, instead of 
yielding FAs, the dark solitons will repel and move away from each other, which is not 
shown here. Hence, coupled dark solitons will transform into FAs if the coupling is in 
an intermediate region. 



We depict in Figs. 2(a) and 2(b) a numerically obtained FAs and its eigenvalues, 
respectively. We observed that the FAs solutions indeed exist stably only for k < 1/3. 
As k approaches 1/3, the amplitude of the humps in the imaginary parts tends to zero 
and the FAs changes into coupled dark solitons. 

Because solitons of ([I]) without a trap are translationally invariant, they can move 
freely in space. This motivates us to study the existence and stability of solitons in a 
moving coordinate frame ^ = x — vt (cf. (E])). 

We present in Fig. 3(a) an FA travelling with v = 0.2 for k = 0.1. One can see 



the deformation in the shape of the soliton due to the nonzero value of v. The shape 
of the deformation suggests the correction ansatz (|8]). For a fixed v, if the coupling 
constant is increased further, the FAs changes into coupled dark solitons at a critical 



value kce < 1/3. The plot of kce as a function of v is shown in Fig. 3(b) as filled circles. 
The solid curve in the same figure is the graph of ([9]), where we obtain perfect agreement. 
We do not show the profile of dark solitons for nonzero v as exact solutions are available. 

After examining the existence of the solitons, next we study their stability. We 
have calculated the stability of FAs for several values of k and find that the FAs are 
always stable in their existence domain. Extending the calculation to coupled dark 
solitons, we found that they are unstable for the values of the coupling constant k < k^e 
corresponding to each value of the velocity v. For the values oik > kce, the dark soliton 
becomes stable. The Runge-Kutta method of order four has been used to verify the 
time dynamics of the FAs as well as dark soliton. 



4-2. Coupled BECs with a trap: Vt ^ Q 

Next, we consider the existence and stability of FAs and coupled dark solitons in the 
presence of a harmonic trap. Two FAs for two different values of coupling constant k are 
shown in Figs. 4(a) and 4(c) with f2 = 0.1. In both panels, we present in dash-dotted 
curves our approximation (ITOl) . In panel (c), we also depict our Gaussian approximation 



f l2T]) . where one can note that the ansatz is slightly better than the Thomas- Fermi ansatz 
f lTOj) . The relatively good agreement deviates rapidly as k increases towards the critical 



coupling. In Figs. 4(b) and 4(d) , we present the eigenvalues of the FAs, where the 
solitons are found to be unstable and this agrees with the result obtained in ( IT9|) . 

The finding that the FAs above are unstable suggests us to investigate whether FAs 
are always unstable for nonzero Vt. The result is summarized in Fig. |5l Shown in solid 
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Figure 4. Numerically obtained FAs solutions of ([TJ and their eigenvalues for = 0.1, 
w = 1, /i = 1, (a-b) k = 0.1, (c-d) k = 0.008. The dasli-dotted curves in (a)-(c) are the 
Thomas- Fermi ansatz (jlOp and the dashed curves in (c) are approximations obtained 
through the variational approach using the Gaussian ansatz . 



curves are the first few lowest squared eigenvalues as a function of k for = 0.1. As 
stability corresponds to > 0, the figure shows that there is a critical value kcs above 
which FAs are stable, which for the parameter value above is approximately 0.145. We 
also present the analytical result (|T9ll in dashed-dotted line, where qualitatively good 
agreement is obtained. 

We also perform the time dynamics of the unstable FAs in Fig. 4(a) Shown in Fig. 
[6] is the typical evolution of unstable FAs, where one can see that the instability causes 
the solitons to oscillate about the minimum of the trap non-sinusoidally. Hence, despite 
the instability, the FAs still persists. This result rather applies generally to other values 
of coupling constants where FAs are unstable. 

As for the existence of FAs in the presence of a magnetic trap, we observe that the 
critical coupling kce above which FAs do not exist is almost independent of the trapping 
parameter Q, which is reasonable as the width of the non-zero imaginary part of the 
FAs solution is small compared to the width of the trapping parameter Q when k — > kce- 

We have also studied the existence and the stability of dark solitons. Shown as 
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Figure 5. The first few lowest squared eigenvalues of FAs (solid) and dark solitons 
(dashed) as a function of fc for fi = 0.1. The dash-dotted curve represents the 
approximation (|19p . Shown is l^pjl- 
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Figure 6. A numerical evolution of the FAs in Fig. 4(a) for fl = 0.1 and k ~ 0.1. 
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Figure 7. Numerical evolution of coupled dark solitons for £7 = 0.1 and k = 0.1. 
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dashed lines in Fig.Oare the first few lowest squared eigenvalues of coupled dark solitons, 
where similarly as before for k less than the critical coupling kce coupled dark solitons 
may transform into an FA as shown in Fig. [71 



5. Conclusion 



We have studied the existence and the stability of FAs and coupled dark solitons 
in linearly coupled Bose-Einstein condensates, both in the absence or presence of a 
harmonic potential. Without a trap, we have obtained travelling FAs and determined the 
range of values of the coupling constant k where they exist. In the presence of a magnetic 
trap, it has been demonstrated that one can obtain a stable FAs by controlling k. In 
the region of k where an FA is unstable, it has been shown that the FAs still exists, but 
oscillates about the minimum of the trap non-sinusoidally. Theoretical approximations 
based on variational formulations have been derived and used to explain the numerical 
results. 

In the future, interaction of multiple FAs and dark solitons in the same setup as 
well as the existence and the stability of higher dimensional excitations, such as vortices, 
will also be investigated and reported elsewhere. 
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